Surrogate Taxa and Fossils as Reliable Proxies of Spatial Biodiversity Patterns in Marine Benthic Communities 


Abstract
Rigorous documentation of spatial heterogeneity (beta-diversity) in present-day and preindustrial ecosystems is required to assess how marine communities respond to environmental and anthropogenic drivers. However, the overwhelming majority of contemporary and paleontological assessments have centered on single higher taxa. To evaluate the validity of single taxa as community surrogates and paleontological proxies, we compared macrobenthic communities and sympatric death assemblages at 51 localities in Onslow Bay (North Carolina, USA). Compositional heterogeneity did not differ significantly across datasets based on live mollusks, live non-mollusks, and all live organisms. Sympatric death assemblages and their subsets were less heterogeneous, likely reflecting time-averaging induced homogenization. Nevertheless, in the context of the scale of variability of beta diversity in other marine ecosystems, all datasets were >80% concordant in all cases. Datasets yielded concordant bathymetric gradients, producing nearly identical ordinations consistently delineating habitats. Congruent estimates from mollusks and non-mollusks suggests that single groups can serve as reliable community proxies. The high spatial fidelity of death assemblages supports the emerging paradigm of Conservation Paleobiology. Integrated ecological and paleontological data employing surrogate taxa can be used to quantify anthropogenic changes in ecosystems and advance our understanding of spatial and temporal aspects of biodiversity.
	
1. Introduction 
Understanding spatial aspects of biodiversity is vital for effective planning of marine protected areas and coastal resource management [1–3]. Consequently, it is essential to identify processes influencing biodiversity on multiple scales [4,5], and document how and why community composition varies within and across habitats (i.e., beta-diversity [?] [6] and related spatially explicit approaches). If ?-diversity is an outcome of environmental conditions, conservation efforts should target the preservation of spatial organization or species-environment relationships, and not resource abundance [7]. In terrestrial systems, ?-diversity has been used widely in conservation and resource management [5,8–12]. However, until recently, few studies have explored ?-diversity in marine communities [13], focusing predominantly on bathymetric diversity gradients in the deep sea. Quantitative studies of ?-diversity in coastal or inner shelf habitats are sparse, and typically examine only a single taxonomic group. Whereas previous studies have provided critical insights into spatial structuring of ecosystems, it remains unclear whether heterogeneities manifested within single taxonomic groups reflect ?-diversity of communities [14]. Yet the use of a single taxon, or biological surrogate, to estimate biodiversity is becoming increasingly important in marine conservation because of the costs of marine biodiversity surveys, undescribed species, and challenging species identifications [15–17]. Although empirical assessments of ‘entire communities’ are not feasible, it is viable to collect multi-taxic data to evaluate the effectiveness of surrogate taxa as proxies for more representative portions of communities.
In addition to spatial data retrievable from present-day ecosystems, the fossil record archives historical data that are often numerically adequate for computing ?-diversity. Consequently, paleontological data derived from the most recent fossil record can potentially augment ecological studies by providing a historical perspective on spatiotemporal dynamics of marine ecosystems over timescales inaccessible to ecology [18–24]. However, the utility of the marine fossil record as a historical proxy of biodiversity is uncertain, as most fossil assemblages are time-averaged mixtures of past communities amalgamated over centennial-to-millennial time scales (e.g., [25–28]). Furthermore, the fossil record primarily preserves organisms possessing biomineralized parts [29–31]. Live-dead studies are widely employed to assess the extent to which dead remains replicate living communities. However, these studies are typically restricted to single groups of organisms. In aquatic systems, fidelity studies have focused primarily on mollusks [18,19,32–35], with only a few case studies dedicated to other groups [36–40]. The emerging consensus points to high live-dead fidelity [35], although ?-diversity may often be somewhat reduced [35] and alpha diversity is typically elevated [32,33,35]. Arguably, despite fossilization filters and time-averaging, mollusk death assemblages capture the underlying variation in molluscan communities. However, mollusks are only a subset of communities. Whether paleontological data are adequate to estimate spatial structuring of communities remains unclear.
Thus, two largely untested assumptions underlie ecological and paleoecological efforts: (1) a single group of marine organisms (typically mollusks) can approximate ecological patterns characterizing communities (implicit to any study that employs a single group of organisms), and (2) high fidelity of one fossil group such as mollusks validates the fossil record as an ecological community proxy. 
To assess these two key assumptions, we examined spatial ecological patterns in macrobenthic marine invertebrate communities from shallow marine environments within Onslow Bay (North Carolina, U.S.A.). Specifically, we estimated heterogeneity in faunal composition across sampled localities using various measures of  ?-diversity: [?VARIANCE] [41,42], [?SHANNON] [43,44], and multivariate dispersion [?DISPERSION] [45,46]. Spatial turnover along the regional water depth gradient was estimated using pairwise comparisons of faunal similarity [?GRADIENT] [42]. To assess fidelity, we compared ?-diversity estimates produced by Live Assemblages [LA], Death Assemblages [DA], Mollusk Live Assemblages [Mollusk LA], Robust Mollusk Live Assemblages [Robust Mollusk LA] that excluded species with low fossilization potential (e.g., shell-less nudibranchs and thin-shelled fragile species such as Anomia simplex), Non-Mollusks Live Assemblage [Non-mollusk LA] with mollusk species removed, and Mollusk Death Assemblages [Mollusk DA]. The comparisons of estimates derived for the entire community [LA] with estimates obtained for DA and increasingly restricted data subsets (Mollusk LA, Robust Mollusk LA, Non-mollusk LA, Mollusk DA) represents a series of direct empirical tests of the fidelity of ecological proxy data provided by the live and dead macrobenthos. Finally, we compare estimates obtained for each dataset against 21 datasets from other ecosystems, to evaluate whether all datasets result in the same comparative outcomes across ecosystems. We employ this comparative approach here to assess the two fundamental assumptions invoked implicitly in many studies of marine communities, both fossil and modern.

2. Methods
We surveyed marine benthic macro-invertebrate communities in coastal and inner shelf habitats near the city of Beaufort, North Carolina, USA (figure S1 and table S1). Samples were obtained by dredging 51 localities over four field seasons (June 2011, November 2011, May of 2012, and April 2013). Although previous studies indicated that seasonal variability in community composition was negligible in this region [47], 31% of our localities were resampled in different seasons to reduce potential seasonal effects and improve locality-level estimates of alpha diversity and relative species abundances. Sympatric DA were collected in conjunction with LA, and macro-invertebrates in DA were identified to the lowest taxonomic level possible and counted (table S1). Samples were pooled to derive locality-level estimates, resulting in two data matrices (LA and DA) of species counts with 51 localities and 179 and 160 species respectively. The restricted matrices were derived from the LA and DA datasets for taxonomic subsets: Mollusk LA, Robust Mollusk LA, Non-mollusk LA, and Mollusk DA (table S2). One locality was removed from all analyses (locality 39), as it did not contain any species for one of the data subsets. The LA included multiple higher taxonomic groups, and adequate representation of infaunal organisms (27% of species were infaunal). Encrusting species were excluded. 
The fidelity of alpha diversity and evenness between the LA and DA was evaluated by measuring LA-DA offsets in evenness/diversity estimates [32,33,35]. Evenness was measured using Hurlbert’s Probability of Interspecific Encounter (PIE), calculated as PIE=[N?((N-1) )](1- ?_(i=1)^S?p_i^2 ), where N=sample size, S=richness, and pi=proportion of species i  [48–50]. PIE’s statistical properties make it particularly suitable for the alpha-level fidelity analyses (51). Evenness was measured as the differences between LA and DA PIE values (? PIE), with ? PIE values of 0 indicating identical evenness estimates for two compared datasets and positive/negative values indicating higher/lower evenness in one the two datasets [33]. PIE [51] and ? PIE [33] are not sensitive to variation in sample sizes or number of samples. Differences in richness were measured as the difference between logarithmic locality-level richness of two compared datasets (? ln S). For each locality richness estimates were rarefied to account for differences in sample size [33] between the compared datasets. 
?VARIANCE, a within-habitat heterogeneity in community composition [42], was estimated for untransformed data as the total community composition variance, calculated here as the total sum of squares of the Bray-Curtis dissimilarity matrix divided by the number of localities less one (n - 1) [41]. ?SHANNON was calculated on untransformed data using Shannon entropy [43,44] and used for direct comparison with estimates from the literature estimates complied in [52]. ?DISPERSION was calculated as multivariate homogeneity of group dispersions [45,46] using scores from Principal Coordinate Analysis of Bray-Curtis similarity matrix obtained for standardized datasets (double-relativization (‘Wisconsin’) standardization of square-root-transformed specimen counts). Spatial median was used to estimate centroids of datasets. All ?-diversity indices were also calculated after sample size standardization. For each dataset 1000 subsampled locality sets were generated using the minimum sample size observed across the six datasets at each locality. The produced distribution provides resampling estimates of statistical uncertainty of the estimated beta diversity measures.
?GRADIENT was measured as pairwise Bray-Curtis similarity between all localities along a depth gradient, to examine change in community structure from one sampling unit to another [42,53]. Water depth was targeted because environmental conditions (e.g., light penetration, disturbance, salinity, temperature, etc.) should vary with water depth and distance from shore, resulting in changes in benthic community composition with depth [54–59], and a strong relationship between community composition and depth for the LA has been demonstrated previously [60]. Exponential decay models were used to evaluate the rate of turnover along the gradient by comparing pairwise Bray-Curtis similarity with pairwise depth difference [42]. Note that pairwise distance here pertains to differences in water depth and not differences in spatial distance. Samples may be as much as 32 km apart, but differ in water depth by <1m (figure S1). However, on average, similarity tends to be high when comparing localities that occur at similar depths, and low when comparing localities with depth differences > 2m, regardless of spatial distance.
To graphically explore changes in similarity along the depth gradient and compare the ordering of samples across different datasets, Nonmetric Multidimensional Scaling (nMDS) was performed using Bray-Curtis similarity for ‘Wisconsin’-standardized, square-root transformed specimen counts. Because nMDS analyses in three dimensions yielded more acceptable stress values ? 0.2 and produced ordinations comparable to nMDS performed in two dimensions, three-dimensional nMDS was used in the final analyses. Environmental variation within- and between-habitats was compared using PERMANOVA [61] of Bray-Curtis similarity values. A total of 9999 permutations (plus the observed value) were performed to generate resampling distribution of pseudo-F statistic to assess the significances of partitions (in this case, habitats). Higher values of F suggest greater multivariate differences between habitats. Significance of differences in ?-diversity between datasets or habitats were tested using permutational (9999 iterations) estimates of homogeneity of multivariate dispersions, which can be thought as multivariate extension of the univariate Leven’s test for homogeneity of variance [45]. Confidence intervals for average distances to habitat centroids expected for habitat-invariant beta disparity was estimated by randomizing localities across habitats (1000 iterations).
All analyses were performed using R (version 2.15.1). Bray-Curtis pairwise distances, ?VARIANCE estimate computations, and iterative sampling were carried out using custom scripts written by the authors. nMDS and ?DISPERSION tests for homogeneity of multivariate dispersions were performed using ‘Vegan’ [62]. ?SHANNON was computed using ‘Vegetarian’ [63].

3. Results 
The LA dataset includes 51 localities, 179 species from 7 phyla (Annelida, Arthropoda, Brachiopoda, Cnidaria, Echinodermata, Mollusca, Porifera), and 11,551 individuals. The DA dataset, sampled at the same 51 localities, totaled 57,116 individuals, 160 species, and the same 7 phyla. Because all other datasets represent subsets of either LA or DA, two types of assessments are possible: (1) non-overlapping comparisons of datasets that represent unique datasets  (Mollusk LA vs. Non-mollusk LA, Mollusk LA vs. Mollusk DA, and LA vs. DA); and (2) overlapping comparisons, where one dataset is a subset of the other (e.g., LA vs. Mollusk LA, or DA vs. Mollusk DA). Non-overlapping comparisons are more appropriate statistically and more conservative when assessing concordances between datasets. However, overlapping comparisons are also relevant because they represent direct analogues of surrogate and fossil datasets (i.e., subsets of unknown larger wholes). Except for the Robust Mollusk LA (discussed only peripherally here), the subset datasets do not vary dramatically in their proportional representations of the entire datasets (table 1): Mollusk LA include 53% of LA species and 33% of LA individuals, Non-mollusk LA include 47% of LA species and 67% of LA individuals, and Mollusk DA include 73% of DA species, and 22% of DA individuals.

3.1. Non-overlapping Comparisons 
For Mollusk LA and Non-mollusk LA datasets (and their sample-standardized variants), all measures of ?-diversity yielded similar estimates (table 1) and the two datasets were statistically indistinguishable in terms of multivariate dispersion (table 3). When compared to 21 literature estimates from other ecosystems (figure 1), the sample-standardized Mollusk LA and Non-mollusk LA datasets were 100% mutually consistent. Specifically, ?SHANNON estimates based on Mollusk LA and Non-mollusk LA were significantly elevated in four case-study comparisons (figure 1 C, K, P, Q) and significantly depressed in 17 comparisons (figure 1 A-B, D-J, L-0, R-U). When assessing spatial structuring of faunal assemblages across habitat, Mollusk LA and Non-mollusk LA datasets produced consistent nMDS ordinations. For both datasets, samples grouped by habitats along a nearshore-offshore depth gradient and convex hulls of habitats were comparable in size and overlap (figure 2 A, B). Also, for both datasets, the habitats were statistically distinct in faunal composition (table 2). Within-habitat estimates of ?-diversity were also mutually concordant, with both datasets indicating that ?-diversity in nearshore habitats was significantly depressed relative to a null randomization model (figure 3). The onshore-offshore sorting of localities (figure 2) points to a depth-related community turnover gradient (?GRADIENT), where similarity in faunal composition of localities is inversely related to difference in water depth between those localities (figure 4). For both Mollusk LA and Non-mollusk LA datasets (figure 4 B-C), the gradients followed comparable trajectories (figure 4 B-C, G) with similar decay rates (figure 4 H), suggesting comparable depth-related ?-diversity gradients. However, depth-invariant offsets in Bray-Curtis similarity values (figure 4 G) and decay curves (figure 4 H) suggest that the Mollusk LA fauna was more homogenous (higher pairwise similarity estimates) than the Non-mollusk LA fauna. Comparisons of coefficients of determinations estimating the relation between depth-difference versus faunal similarity for all datasets indicate that consistent proportions of variation are explained by the onshore-offshore gradient (figure S2). Finally, average alpha diversity and evenness of localities were on average slightly higher for Mollusk fauna (figure S3 A), but the difference was not significant (table S4). 
The LA vs. DA comparison was, in many aspects, similar to the Mollusk LA vs. Non-mollusk LA comparison. Pairwise comparisons of the DA-based ?SHANNON estimates against literature estimates were mostly consistent with LA-based outcomes (figure 1; 81% concordance rate). LA-based and DA-based multivariate ordinations were concordant (figure 2 C, D) and mutually consistent in statistical outcomes (table 3). For both datasets, within-habitat estimates of ?-diversity indicated significantly depressed ?-diversity in the nearshore habitat (figure 3). However, in contrast to the Non-mollusk LA – Mollusk LA comparison, ?-diversity of DA was significantly lower comparing to LA (table 3). Also, DA-based ?GRADIENT estimates were offset more substantially relative to LA (figure 4 G), and the DA decay curve was less steep than the LA curve (figure 4 H), suggesting a stronger spatial homogenization within and across habitats. Finally, DA-based alpha diversity and evenness estimates were significantly elevated compared to LA (figure S3 B, table S6).
Mollusk LA and Mollusk DA estimates of beta diversity were mutually concordant, including statistically indistinguishable estimates of ?DISPERSION (table 3), 100% concordance in comparisons with literature case studies (figure 1), comparable ordinations (figure S5), and consistent indication of significantly depressed ?-diversity in the nearshore habitat (figure 3). One exception are harbor samples, for which Mollusk DA values were significantly depressed relative to Mollusk LA estimates. However, this offset is likely an artifact due to the very small number of harbor localities. Also, Mollusk LA and Mollusk DA produced comparable, if slightly offset, ?GRADIENT trajectories (figure 4 G), and consistent decay rates (Figure 4 H). Thus, Mollusk DA and Mollusk LA were mutually more congruent in terms of ?-diversity than was observed in the LA-DA comparisons. However, as in the case of LA vs. DA, alpha diversity and evenness were significantly elevated for Mollusk DA compared with Mollusk LA (figure S3 table S6). 
	The largely concordant patterns between the Non-mollusk LA and Mollusk DA (figures 1-4 and S4, tables 1-2 and S6) demonstrate that compositional fidelity is not a prerequisite for spatial fidelity despite no shared species between these two assemblages. Compositional fidelity therefore should not be considered a prerequisite for spatial fidelity (an in-depth discussion of fidelity will be presented elsewhere).
3.2. Overlapping Comparisons
The LA dataset is the most comprehensive estimate of the spatial structuring of local ecosystems of Onslow Bay, and thus, can serve as a comparative baseline for assessing the performance of surrogate taxa (Mollusk LA, Non-mollusk LA, Robust Mollusk LA) and paleontological proxies (DA and Mollusk DA). Except for Robust Mollusk LA (which represent a small fraction of all LA data), the two subsets provided reasonable approximations of all aspects of ?-diversity and spatial structuring of faunal assemblages across and within habitats (figures 1-4; tables 1), although Mollusk LA displayed slightly more depressed estimates of ?-diversity than Non-mollusk LA (figure S6) and, consistently, stronger offsets in the ?GRADIENT (figure 4 G, H).
Paleontological proxies of LA, produced ?-diversity estimates that were more strongly suppressed (figure S6; figure 4 G, H) and significantly overestimated average alpha diversity and evenness of individual localities (figure S3). However, comparing to DA, Mollusk DA produced estimates of ?-diversity that were more consistent with LA-based estimates (figures 1 and 4; figure S6; table 3), including a better agreement with literature comparisons (figure 1; Mollusk DA 90.4% concordant, DA 76.2% concordant) and ?GRADIENT estimates (figure 4 G, H). Nevertheless, both DA and Mollusk DA produced multivariate ordinations consistent with LA (figure 2, figure S5) and were qualitatively congruent in many other aspects of the analysis (figures 3 and 4; table 2). 

Discussion 
The high concordance of ?-diversity metrics and nMDS ordination plots observed when comparing non-overlapping (Mollusk LA versus Non-mollusk LA) and overlapping (Mollusk LA and Mollusk DA versus LA) datasets indicates that spatial structuring of the macrobenthic community in the Onslow Bay area can be reliably approximated when data are restricted to one group of benthic organisms. Different faunal components of the community are structured in a spatially consistent manner. However, although habitat segregation is significant for all datasets, the LA, and the Non-mollusk LA analyses produced a slightly weaker separation of the habitats (LA F=2.94; Non-mollusk LA F=2.77) than the Mollusk LA dataset (F=3.31). That is, in the case of Onslow Bay, mollusk-based data differentiate habitats more clearly than non-mollusks or the entire LA (this difference cannot be attributed to differences in taxonomic resolution, given that both mollusks and non-mollusks include comparably small fractions of unresolved taxa:23% and 11%, respectively).  In contrast, Non-mollusk LA performed subtly better in reproducing LA-based beta diversity gradients and overall levels of spatial heterogeneity in faunal composition (figure 4). However, the weaker correlation between the Non-mollusk LA and LA may indicate much higher within-habitat discordance (figure S4), possibly a scaling artifact. Finally, any biases introduced by the use of surrogate taxa to assess spatial structuring of a given ecosystem are meaningful only in cases when they produce inconsistent outcomes in comparisons to other ecosystems. Mollusk LA and Non-mollusk LA datasets were mutually consistent in 100% of comparisons to literature estimates from other ecosystems, and both those surrogate datasets were also 90% concordant with LA dataset (figure 1). Only in 1 out of 21 comparisons, the discrepancies between the LA and surrogate datasets mattered. In summary, although patterns in ?-diversity vary among ecosystems and organisms [5,64], the outcome of this study is consistent with previous studies demonstrating that subset-taxa surrogates perform well [15], and that cross-taxa surrogates provide congruent representations of compositional turnover [16]. These results tentatively indicate that spatial ecological analyses restricted to one group of organisms may produce meaningful proxies for assessing spatial patterns estimated by multiple groups, a potentially powerful tool in conservation planning [17]. 
The two paleontological proxy datasets (DA and Mollusk DA) also produced ordination plots that are in agreement with the corresponding estimates obtained for live communities. Thus, despite post-mortem transport, time-averaging, and taphonomic filtering of species with low fossilization potential, compositional differences between habitats observed in living communities are reflected in the Onslow Bay death assemblages. The suppressed estimates of ?-diversity observed in DA are consistent with the hypothesis that death assemblages may be affected by spatial homogenization because of time-averaging and post-mortem transport [28,35]. The slightly elevated locality-level diversity in DA assemblages is also consistent with this explanation: spatial and temporal mixing is not only expected to suppress ?-diversity but often results in inflated alpha diversity (table 4 in [35]). However, the bias in ?-diversity estimates observed for paleontological proxy datasets is largely trivial when considered in the context of the magnitude of variation in ?-diversity estimates observed across ecosystems (figure 1). Compared to the entire life assemblage (LA), both DA and mollusk DA produce consistent conclusions in relative comparisons to other ecosystems (>80% concordance in both cases). In general, no matter which of the live or dead datasets is used, the same conclusion is reached: the spatial heterogeneity in the Onslow Bay area is low comparing to previously reported estimates from other ecosystems. The results are thus consistent with the growing evidence for the high informative value of death assemblage data [28]. This analysis suggests that death assemblages are an adequate proxy for the spatial structuring of the entire community and not just its mollusk component. The fact that Mollusk DA datasets produced less biased estimates than the entire DA is particularly encouraging given that mollusks are the most common target of conservation paleobiology and fidelity studies. It is important to note, however, that mollusk assemblages tend to have relatively high live-dead agreement in species identity and rank-order abundance [19,34,65], and as a result, may record spatial structuring of ecosystems more accurately than other groups (e.g., arthropods, echinoderms) that have not been examined adequately so far and are insufficiently represented in our data to allow for such analyses.
Three alternative explanations can be proposed for the relatively high spatial fidelity of death assemblages. First, the high dead-live concordance may indicate that death assemblages reliably record living communities. If this explanation is correct, then the spatial structuring of local macrobenthic communities have not changed dramatically despite prolonged anthropogenic activities that have affected Onslow Bay [66–70], and the regional ecosystem structure has not shifted notably from its historical state. Second, the high fidelity may indicate that the regional death assemblage has already shifted to a new state and also records altered ecosystems. This is unlikely given that multi-centennial to multi-millennial time-averaging is the norm for surficial shell accumulations in comparable, sediment-starved coastal and shelf settings [71,72], as demonstrated by direct dating of surficial shell remains in other regions [27,28]. Finally, changes induced by human activities may have led to community homogenization, a common outcome of anthropogenic changes [73], which would make living communities more comparable to time-averaged records archived in death assemblages. However, the observed compositional differentiation of habitats (figure 2) and the beta diversity gradient (figure 4) indicate consistently that the regional system has not been homogenized. The results are most congruent with the first explanation: death assemblages provide qualitatively meaningful archives of the spatial structuring of living communities.
An effective mitigation of the current marine biodiversity crisis, attributable to widespread human activities in marine habitats [74], is not viable unless we develop a robust understanding of both spatial structuring of present-day ecosystems [5,52] as well as historical dynamics of past communities, prior to, and after, the onset of anthropogenic changes [18,21]. The positive outcomes of this study reinforce the validity of spatial ecological estimates based on surrogate taxa. The results also indicate that paleoecological approaches can provide adequate historical measures of ?-diversity. By integrating ecological and paleontological data for surrogate-taxa, we can advance our understanding of spatial and temporal patterns in biodiversity, investigate the long-term dynamics of past communities, and quantify the magnitude of recent changes in spatial structuring of marine ecosystems.